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We present a Monte Carlo wavefunction method for semiclassically modeling spin-| systems in 
a magnetic field gradient in one dimension. Our model resolves the conflict of determining what 
classical force an atom should be subjected to when it is in an arbitrary superposition of internal 
states. Spatial degrees of freedom are considered to be an environment, entanglement with which 
decoheres the internal states. Atoms follow classical trajectories through space, punctuated by 
probabilistic jumps between spin states. 

We modify the conventional Monte Carlo wavefunction method to jump between states when 
population transfer occurs, rather than when population is later discarded via exponential decay. 

This results in a spinor wavefunction that is continuous in time, and allows us to model the classical 
particle trajectories (evolution of the environment variables) more accurately. The model is not 
computationally demanding and it agrees well with simulations of the full spatial wavefunction of 
an atom. 


INTRODUCTION 

A semiclassical simulation is one in which some degrees 
of freedom of a system are treated quantum mechanically 
and some classically. 

Semiclassical simulations in which the positions 
and momenta of atoms are treated classically enjoy 
widespread use in cold atom physics, for example in laser 
cooling m and magnetic trapping/evaporative cool¬ 
ing 0- The benefit of not simulating the positions and 
momenta quantum mechanically is obvious: one need not 
spend the considerable computational resources required 
to do so. In addition, methods that define explicit classi¬ 
cal trajectories are readily amenable to simulating colli¬ 
sions amongst ensembles of spins in magnetic fields using 
molecular dynamical methods such as direct simulation 
Monte Carlo (DSMC) G . 

The downside is that superpositions of (and hence in¬ 
terference between) position and momentum states can¬ 
not be included in such a simulation. This is usually not 
a problem, as the position evolution of an atom in a laser 
cooling simulation for example is well approximated by 
the action of a classical force—which doesn’t result in 
superpositions of position states. 

Furthermore, although atomic wavepackets tend to ex¬ 
pand over tirncQ incoherent atomic scattering results in 
frequent position measurements such that atoms in a 
thermal cloud can be described by Gaussian wavepack¬ 
ets of a well defined average size that depends on the 
collective properties of the cloud [8.. 

So it would seem that such semiclassical simulations 
are ideal for cold atom simulations. When are they not? 


1 For a thermal gas, they expand at a rate equal to the thermal 
velocity [6], which for a room temperature rubidium atom results 
in a wavepacket the size of a grapefruit after one millisecond! 


Semiclassical simulations run into problems when it isn’t 
apparent what classical force to use. If an atom is in a 
superposition of two states with different magnetic mo¬ 
ments in a magnetic field gradient, what force should it 
be subjected to? Corresponding to which magnetic mo¬ 
ment? The answer is of course both, and an actual atom 
in such a situation will evolve into a spin position entan¬ 
gled state. 

Since we wish to continue using semiclassical meth¬ 
ods, what can we do about this? One approach is to 
use the average force that the atom should feel, and this 
will indeed result in the correct expectation value of po¬ 
sition over time (Ehrenfest’s theorem 7]). This works 
well when the atoms are mostly in one eigenstate, for ex¬ 
ample in Sisyphus cooling 01 where superpositions are 
dominated by groundstate population due to far detuned 
lasers and are short lived due to spontaneous emission 
(which puts atoms into eigenstates once more). 

This approach does not work however when there 
are non-negligible and long-lived superpositions between 
states that see different potentials, which brings us to our 
motivation for developing the method that is the subject 
of this paper: forced evaporative cooling of neutral atoms 
in a magnetic trap. 

In a magnetic field, the internal energy eigenstates of 
an atom are the states with well-defined spin projection 
in the local direction of the magnetic field, n = jgy. 
Atoms that are spin-aligned with the local field are at¬ 
tracted to regions of low field strength, and atoms anti¬ 
aligned are attracted to high field strength (or the reverse 
if the atom has a negative Lande ^-factor). As such, 
atoms can be trapped in magnetic field configurations 
that have a spatial minimum of magnetic field strength 
somewhere, such as the quadrupole trap 0 often used in 
Bose-Einstein condensation mm- 

However, if the magnetic held direction is varying in 
space, these locally spin aligned or anti aligned states 
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are not eigenstates of the Hamiltonian an atom sees as 
it moves. If the atoms are moving slowly or the rate of 
change of the field direction is small, then the spins adia- 
batically follow the field and remain in a local eigenstate. 
However if they move more quickly, or though a region of 
rapidly changing field direction, transitions between local 
eigenstates can occur, called Majorana transitions ill^j . 
This leads to atoms that were previously trapped be¬ 
coming untrapped, and to losses in magnetic traps. Such 
Majorana losses are a problem in cold atom experiments, 
as cold atoms spend more time near the center of the trap 
where field directions change rapidly. This leads to both 
atom losses and overall heating of the atom cloud. 

Using Ehrenfest’s theorem leads to qualitatively incor¬ 
rect results in this situation. An atom, having passed 
near to the field zero and transitioned into say, a 60:40 
superposition of trapped and anti-trapped states, will feel 
a weakly trapping force according to Ehrenfest’s theorem. 
It will not lose as much kinetic energy departing the field 
zero as it gained approaching it, and will eventually pass 
this kinetic energy on to other atoms in collisions. 

An actual atom in this situation will of course diverge 
into two trajectories. The component of the superposi¬ 
tion in the trapped spin state will remain on as tight an 
orbit about the field zero as before, not gaining any net 
energy from its close encounter. The anti-trapped com¬ 
ponent, conversely, will accelerate away from the field 
zero, and at typical evaporative cooling densities not col¬ 
lide with any other atoms on the way out. 

The archetypal example of spin-position entanglement 
is the Stern-Gerlach experiment |13) . in which spin-| 
atoms—initially in an equal superposition of two spin 
projection states—traverse a magnetic field gradient. 
Two distinct trajectories are observed resulting from the 
two distinct forces on the two spin projection states. 
The expectation value of the force, however, is zero, and 
so simulating the Stern-Gerlach experiment with Ehren¬ 
fest’s theorem similarly results in the qualitatively incor¬ 
rect result of a single trajectory, in the center of where 
the two distinct trajectories would be. 

We solve this problem using a method in which there 
is only one classical force (corresponding to a single 
spin projection state) on each atom at each moment in 
time, despite the possibility of arbitrary superpositions 
of internal states. The method can be visualized by 
considering two spin wavepackets separating as per the 
spatial Schrodinger equation—with the spin components 
accelerating in opposite directions—and asking “What 
spin-state populations do I see if I follow the spin up 
wavepacket?”. As shown in figure [T] one sees the spin 
down population decreasing over time until only the spin 
up population remains. The rate that it does so depends 
on how fast the wavepackets separate. Likewise, if one 
follows the spin down component, one sees the spin up 
population decrease to zero. Our method does just this, 
following the trajectory of one spin component and de¬ 


caying the other. 

Although atoms mid-spin-flip during evaporative cool¬ 
ing do indeed evolve into spin-position entangled states, 
once separated, we are nonetheless happy to neglect pos¬ 
sible future interference between these states. Firstly, 
the probability of the two wavepackets coming together 
again is low, and secondly, collisions with other atoms 
will soon occur, entangling atomic positions and decreas¬ 
ing the probability of future wavefunction overlap expo¬ 
nentially in the number of collisions. Lastly, even if spa¬ 
tial interference does occur, it has little effect on our sim¬ 
ulation. Interference fringes will be apparent on length 
scales close to the thermal wavelength, on which scale 
moving atoms around slightly won’t affect the collective 
properties of the cloud. For these reasons, classical posi¬ 
tions and momenta are suitable. 

In the next section, we detail how our model is im¬ 
plemented, and compare it to more conventional MCWF 
methods. 


SPIN DECOHERENCE FROM POSITION 
SEPARATION 

Rather than computing the entire density matrix at 
every point in time, Monte-Carlo wavefunction (MCWF) 
methods HM3 track a pure quantum state represented 
by a wavefunction. This is computationally simpler, and 
results in the same statistical outcomes as the full density 
matrix approach. It also has the appeal that pure quan¬ 
tum states are what we actually get when we perform 
quantum measurements, so every run of a MCWF sim¬ 
ulation corresponds to an actual possible experimental 
outcome (20) . 

Decoherence f2TU2TT| is introduced in the MCWF 
method by decaying some components of the wavefunc¬ 
tion and not others HU Eg. Which basis this is done 
in depends on the system-environment interaction, and 
which states to decay and which to keep is probabilistic. 
This arises from considering certain environmental states 
to be ‘classical’, in the sense that one can ignore interfer¬ 
ence effects between them. This is usually valid due to 
the large number of degrees of freedom in environments. 
For example, the MCWF method is most often applied 
to decay of excited atomic states via spontaneous pho¬ 
ton emission |25j . The number of photons in the emitting 
mode is considered to be well defined at all times, as if a 
projective measurement of the photon number were being 
repeatedly performed fl8] . 

Spin decoherence in our system arises from the fact 
that different spin states are subject to different spatial 
potentials, and hence spatially separate over time. The 
spin degree of freedom becomes entangled with the posi¬ 
tional degree of freedom, and future interference between 
different spin states becomes unlikely. Our model ap¬ 
proximates this process as irreversible, which allows us 
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FIG. 1. (Top) Single classical trajectory (dashed line) and 
spinor wavepacket extent (red and blue shaded). (Bottom) 
Spin populations (projected on local B (z)) along classical tra¬ 
jectory. This schematic not to scale. Five relevant points in 
time are labeled, to- The atom is in a 50:50 superposition 
of locally spin up and spin down states. The spin up and 
spin down wavepackets accelerate away from each other. t\\ 
At the center of the spin down wavepacket, we now see a re¬ 
duced spin up population, due to the two spin components 
separating. £ 2 : After sufficient separation time, we see no 
spin up population. 1 3 : The magnetic field vector changes di¬ 
rection rapidly, putting the atom into a 70:30 superposition of 
spin up and spin down. Our model chooses with 70% proba¬ 
bility whether or not to start following the classical trajectory 
corresponding to the center of a spin up wavepacket. In this 
example it chooses to do so. £ 4 : The two wavepackets once 
again separate, and from the center of the spin up wavepacket, 
we see a reduced spin down population. £ 5 : After sufficient 
time we see only spin up population remaining. Note: In 
this schematic we’ve shown the decay in wavepacket overlap 
similarly to how it would actually look in reality: decreas¬ 
ing slowly at first, and then more rapidly as the wavepackets 
accelerate away from each other. However in our model we 
approximate simple exponential decay, and so the curves be¬ 
tween to and £ 1 , and between £3 and £4 are instead exponential 
in shape. 

to model it similarly to environmentally induced deco¬ 
herence. This spatial separation results in a decreasing 
overlap between the wavepackets of the two spin compo¬ 
nents with time, decreasing the contrast of future inter¬ 
ference between the spin populations, eventually to zero 
when they are completely separated. 

This decrease in wavepacket overlap does not hap¬ 
pen at a constant exponential rate, due to the Gaussian 
wavepacket shape and the fact that spin components ac¬ 
celerate away from each other. In order to approximate 
the resulting decoherence as exponential decay, we com¬ 
pute an average time r taken for wavepackets to sepa¬ 
rate, and use its reciprocal as the decoherence rate. This 
is the Markovian approximation often made in decoher¬ 
ence models |25] I2S] ■ It comprises the assumption that 


the environment has no memory; in our case that our 
simulation will not keep a record of how long wavepack¬ 
ets have been accelerating apart. The calculation of r is 
detailed below. 

In the language of decoherence, we have a system to 
be modeled quantum mechanically—the spin state of an 
atom, and an environment to be modeled classically— 
the position of said atom. Our ‘classical’ environmental 
states will be minimum-uncertainty Gaussian wavepack¬ 
ets [Bj. For a spin F atom there are 2 F + 1 such states— 
one for each spin projection eigenstate onto the local 
magnetic field. 

Any interaction between a system and an environment 
will tend to diagonalize the system’s reduced density ma¬ 
trix in the eigenbasis of the system-environment interac¬ 
tion Hamiltonian l22j . which for an atom in a magnetic 
field is the Zeeman Hamiltonian: 


= —4 • B(r), 

(1) 

with eigenvalues: 


E m n = B(r) , 

(2) 


where m n is the spin projection quantum number in the 
direction of B, f is the position operator, gp is the atom’s 
Lande g-factor for the states with total spin quantum 
number A, and /xb is the Bohr magneton. The eigen¬ 
states of this Hamiltonian are spin-position entangled 
states: it acts to separate spin states spatially so that 
they cannot interfere, and it will do so in the spin basis 
of the local magnetic field. As such, this is the basis in 
which we will be decohering states with respect to each 
other. 

It might seem odd to consider the position degree of 
freedom of an atom to be an environment, capable of 
performing measurements on a system. However for an 
interaction to appear as measurement-induced decoher¬ 
ence, it need only cause entanglement that makes future 
interference between the states that it measures negligi¬ 
ble. Quantum erasure experiments for example [27] do 
not illustrate such interactions because the entanglement 
is reversible: interference can be restored. Entanglement 
that is near-complete and irreversible however, can be 
rightly regarded as strong measurement [221 . 

For full, large scale environmental decoherence, such 
interference is unobservable due to the sheer number of 
degrees of freedom in large environments; the probability 
of two fragments of amplitude sharing a common origin 
ever encountering each other again in Hilbert space is 
vanishingly small [22]. 

For us the probability is still small, but not nearly 
as vanishingly so. However we are happy to make the 
approximation: that two of our Gaussian wavepackets, 
once separated, will not coincide in phase space at any 
time in the future prior to a collision (at which point full 
large scale decoherence with the rest of the cloud makes 
our approximation near exact). 
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Take a spin-U system in a magnetic field gradient, for 
which spin and position are entangled. Its wavefunction 
is given by: 


F 

i*>= E | tpm n ) 5S |m„), (3) 

m n ——F 

where m n is the spin projection quantum number along 
the local magnetic field, and (r|^> mn ) = 7/> mn (r) is the 
spatial wavefunction for that spin state (taken to be 
a Gaussian wavepacket following a classical trajectory). 
The total density matrix for a spin-4 system is therefore: 


Gaussian wavepackets. This is not the simple exponen¬ 
tial decay required by the Markovian approximation. We 
seek an average rate of separation t.- 13 instead, such that 
the overlap decreases as: 


= exp 



(9) 


Using ^ as the probability that the wavepackets are still 
overlapped, we may compute the expectation value of the 
time they are no longer overlapped. This is what we will 
use as the time constant t,; 7 in ®: 


P ( r ) 


|c t | 2 |V’tl 2 

|c4,| 2 1-04-I 2 


(4) 


where the il> mn (r) are written as ij)f and for brevity. 
This is the density matrix for the entire system plus en¬ 
vironment we are considering, and so is a pure density 
matrix. If we now do a partial trace over the environmen¬ 
tal degrees of freedom [23], we get the (possibly mixed) 
reduced density matrix for the atom’s internal state: 


Pspin — 


p{ r) dr = 


|c t | 2 c T c|(^ t |V>4.) 

c t c t ('04-IV’t) l c f-| 2 


• (5) 


i,\ _ /o°° i d t 

Tlj {)ij f 0 °°mmmdf 

This gives us: 


7~ij exact — 


y/Zjrhexp [?y 2 ] erfc (V2rj) 
maaijKi (rf) 


( 10 ) 


( 11 ) 


where erfc is the complementary error function, I\ i 
is a modified Bessel function of the second kind ancl 
?7 = m 2 a 3 aijh~ 2 is a dimensionless parameter. 

This is quite a mouthful, and is numerically difficult to 
evaluate. However, in the small wavepacket limit 77 1, 

in which position separation dominates (|8]), one obtains: 


Now consider an atom that starts at t = 0 in a su¬ 
perposition of spin states, but with no spin-position 
entanglement- all the spatial wavefunctions V , m„ (which 
are really functions of time too, as we will soon see) are 
equal. So each of the inner products above is equal to 
unity and we have a pure density matrix for our spin 
degree of freedom. 

As time elapses, the spatial wavefunctions, each corre¬ 
sponding to a different spin state, begin to separate, mov¬ 
ing in different directions. The center of each wavepacket 
will move with the classical force: 

Fm n = -V(-p-B) (6) 

=-g F m n ii B V |B|. (7) 


We model these spatial wavepackets as Gaussians with a 
fixed width. The wavepacket overlap decreases with time 
as they spatially separate^ 


KV>i(i)l ipj(t))\ = exp 


1 

32cr 2 


a 2 / 


2 2 

m a 9 9 

- 9~ a iit 

2 h 2 13 


( 8 ) 


where m is the atom’s mass, aij is the magnitude of 
the relative acceleration between the two wavepackets as 
computed from the force 0 , and a is the size of the 


2 See supplementary material for details. 


(27T 2 )4 I a 

= 2T(|) \/«b 


( 12 ) 


« 1.163 .pF. (13) 

V 

In the large wavepacket limit 77 1, in which velocity 

separation dominates, one instead obtains: 


7~ij vel — 



and the following expression: 


1 


1 


1 


1 

3 


7~ij approx 


T 3 

IJ pos 


T 3 
ij vel 


(14) 


(15) 


gives a good approximation to for all wavepacket 
sizes <7. We have studied the impact of using different 
values of r and found that our method is relatively in¬ 
sensitive to the precise value calculated here. 

We now have a set of decoherence times Tij between 
each pair of states. Note that these are functions of 
the local magnetic field gradient and temperature—they 
should be recalculated appropriately for each atom as the 
simulation proceeds. 

These times describe how long a spin superposition can 
live before the two spin states spatially separate and thus 
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can no longer interfere. As such, in a short time At, our 
pure reduced density matrix will become slightly mixed: 


Pspin( At) 


|c t | 2 c T cJe 

Lcjc^e-^ Icj.! 2 


(16) 


where r = T|j, = (for a spin-1 system there is only one 
spin decoherence time). In this way one obtains a reduced 
density matrix with decaying off diagonals, telling us that 
this is the basis in which spin decoherence occurs, and at 
what rate it does so. 

Conditioned on a particular outcome of a spin mea¬ 
surement, we can then make statements about what the 
wavefunction looks like over time. If the particle is even¬ 
tually determined by measurement to be spin up, then 
the wavefunction over time will have the spin down com¬ 
ponent decaying to zero. This is equivalent to decompos¬ 
ing the reduced density matrix at each moment in time 
into the sum of pure density matrices: 


Pspin(At) 


- |c T | 2 c T cJe-^ - 

_c lC *e-¥ e- 2 ^|cq| 2 _ 

ro 0 I 


+ 


0 (1 — e -2 ^)|cj.| 


(17) 


and discarding the second term each time. The conven¬ 
tional MCWF method allows for the possibility of in¬ 
stead discarding the first term and keeping the second 
one, resulting in an instantaneous change to the pure 
wavefunction being simulated; a ‘quantum jump’. At 
each timestep it chooses between the two above terms, 
weighted by the squared amplitudes of their correspond¬ 
ing wavefunctions. In our method, however, we have al¬ 
ready decided in advance that we will take the first term, 
based on earlier population transfer between the two spin 
states. Alternately, a simulation with our method may 
begin, based on population transfer during some inte¬ 
gration timestep, to start ‘tracking’ the spin down state, 
after which the decomposition is instead: 


Pspin(At) — 


+ 


e 2 "|c t | 2 C|c|e 

_ c^eT^ |cj.| 2 _ 

'(l^e- 2 ^)| Ct | 2 O' 

7 

0 0 


(18) 


once again discarding the second term at each timestep. 
Our method of deciding in advance which term to dis¬ 
card, described in the next section, results in correct 
statistical outcomes for spin measurements (as shown 


in the results section by comparison with the spatial 
Schrodinger equation), and is in that regard identical to 
the conventional MCWF method. However, by bring¬ 
ing forward the decision as to which term of the density 
matrix decomposition to discard, we avoid instantaneous 
changes in the wavefunction, and more accurately simu¬ 
late the classical trajectory by acknowledging a change 
in the spin state sooner. 


METHOD 


Here we present our model in one dimension and for a 
spin-half atom in one dimension. 

In our model, each atom’s state is fully described by 
its position z, its velocity v z , its locally spin up and spin 
down populations c-j- and cp and its currently tracked 
state, that is, which of the internal states |f) or |f) we are 
using to compute the classical force on the atom at that 
moment in time. We use language like ‘we are tracking 
the |f) component’. The currently tracked state can be 
stored simply as an integer, for example +1 for locally 
spin up, and —1 for locally spin down. A spin flip, in the 
context of this method, is when we change, in the course 
of a simulation, which of the two states we are tracking. 

Below, whenever we mention without qualification spin 
up and spin down, |f) and |f), or their amplitudes oj- 
and Cj,; we are referring to states of well defined spin 
projection in the direction of the local magnetic Held the 
atom in question sees, not in the z direction or any other 
lab basis. 

The model is based on tracking one spin component at 
a time, and assuming its spatial wavefunction is a fixed- 
size Gaussian wavepacket with mean position following a 
classical trajectory according to the potential that spin 
component experiences. As in the conventional MCWF 
method, we exponentially decay the other, untracked, 
spin state according to the decoherence time between the 
pair of states. 

Unlike the conventional MCWF method, we do not 
wait until the simulation discards population via expo¬ 
nential decay to consider ‘jumping’ to other states. In 
developing the present model, we observed that doing so 
results in correct spin flip probabilities, but slightly incor¬ 
rect classical trajectories. We attribute this to the time 
delay between population transfer taking place and a 
subsequent ‘jump’ occurring in the conventional MCWF 
method, a delay of on average one exponential time con¬ 
stant. The atom experiences the wrong classical force in 
the meantime, which leads to discrepancies between the 
resulting trajectories and those expected from full spatial 
Schrodinger equation simulations. 

Instead, we determine spin flip probabilities based 
on population transfer between spin states during each 
timestep. As a result, we do not need to ‘jump’ to an¬ 
other state by instantaneously changing the wavefunction 
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when we decide to start tracking a different spin state. 
We simply change which state the simulation is expo¬ 
nentially decaying. This represents a significant depar¬ 
ture from the jumps in the conventional MCWF method, 
the physicality and interpretation of which is the topic of 
more general ongoing debate and discussion mm- 

The instantaneous jumps of the conventional MCWF 
method may be interpreted in the context of our method 
as accounting for antecedent population transfer (on av¬ 
erage, one decoherence time ago). While such simulations 
produce correct final results for the system described by 
the wavefunction, they do not allow correct simulation of 
the corresponding classical environmental states—which 
we are interested in too. To clarify, we have not done 
away entirely with jumps in the non-environment vari¬ 
ables, this would clearly be unphysical since such jumps 
are very real in many experiments [251129] . Our system 
wavefunction still ‘jumps’ from one spin state to another, 
but it does so continuously over the decoherence time, 
instead of instantaneously. For short decoherence times, 
this nonetheless may appear as a practically instanta¬ 
neous change in the system’s wavefunction. 

From the vector spin operator S = { S X ,S V ,S Z }, we 
can construct a spin projection operator S n = n • S in a 
direction described by the unit vector n = {n x , n y , n z }, 
which for us is the direction of the local magnetic field at 
each atom’s position. For a spin-^ system, the projection 
operators onto the eigenstates of S n are: 


and: 


Pt = It) (tl 


= | m n = + §) (m n = + |], 

\(l + n z ) \(n x -iriy) 


Pt = 


L 2 


2 ( ^ 


T in^) 2 (1 n-z) 


Pi = It) (tl 


Pi. = 


= | m n = (m„ = - \\, 

2 (f ^z) in v 

2 ~b ifiy ) 2 (f T n z ) 


(19) 


( 20 ) 


where P^ and P^ are the matrix representations of P^ and 
Pi in the z spin basis. These projection operators are 
used at various steps in the model to project the spinor 
wavefunction onto the eigenstates of spin projection onto 
the local magnetic field. 


conditions intended to approximate those of a spatially 
extended wavepacket, independently draw positions and 
velocities randomly from the position and velocity prob¬ 
ability distributions of the spatial wavefunction in ques- 
tioijfj 

To simulate any initial spin superposition, set c-j- and 
Cj. for each atom accordingly, ensuring normalization is 
satisfied: |c-f| 2 + |cj.| 2 = 1. Compute the eigenvectors of 
the local spin operator S n ; |f) and |j.) in the z basis for 
each atom (or whichever basis you will be simulating in), 
and construct that atom’s initial spinor wavefunction: 

lx) = c t |t) + c i It) • (21) 

For each atom, choose which state initially to track. 
Do so by choosing randomly from the locally spin up 
and spin down states, weighted by their populations |c|| 2 
and |cj.| 2 . For example, if all N atoms start in a 50:50 
superposition of spin up and down, this should result in 
( N± VN )/2 atoms initially being tracked as spin up, and 
the rest as spin down. 

The below steps describe integration for one atom in 
the model. Atoms evolve independently, and so the 
method is trivially parallelizable. 

1. At the start of each integration step, note the state 
populations: 

n t (t) = \c t (t)\ 2 = (x(t)\Pt\x(t)), ( 22 ) 

and: 

n l (t) = \c i (t)\ 2 = (x(t)\Pi\x(t)). (23) 


2. Do ordinary Hamiltonian evolution of the spin state 
for one timestep At: 

I X(t + At)) = exp |x(t)), (24) 

where H is the Zeeman Hamiltonian — p.-B(,z). Any 
integration method can be used, this unitary, first 
order method is shown as an example. 

3. Evolve the position and velocity for one timestep 
At according to the classical force on the state be¬ 
ing tracked: 

F m ntracked = —'V(— \i-B(z)) 

= C/FTTln trackedMB^ l-B^) | • 

Depending on the integration method, this might 
be done simultaneously with step [2] 


Algorithm 

To simulate classical initial conditions, simply set the 
particles’ positions and velocities accordingly. For initial 


3 If you require position-velocity correlations for your initial condi¬ 
tions, use an appropriate joint probability distribution for z and 
v z , the correct derivation of which is beyond the scope of this 
paper. 
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4. Note the state populations again, and determine 
from the state populations noted in step |T] how 
much probability moved from the state currently 
being tracked, to the other state during this inte¬ 
gration timestep. Compute this as a fraction of the 
currently tracked state’s population. For example, 
if the simulation is currently tracking the atom’s 
spin up state, then we have: 


where P U ntracked , either P t or P±, defined in (|19j) 
and (201 is the projection operator onto the state 
not being tracked, and r, given in (15) is the spin 
decoherence time. Again, the numerical method 
used for this integration step is not important, a 
simple Euler method step is shown here as an ex¬ 
ample, and was accurate enough for our examples, 
for which At << r. 


Pa ip 


K(f + Af)| 2 -| Ci (f)| 2 
|c t (f + At) | 2 


(26) 


5. Draw a random number between zero and one. If 
it’s less than pm p , then the atom is to undergo a 
spin flip, if doing so is energetically allowed. The 
above quantity pfn p might be less than zero, in¬ 
dicating probability flow the other way (from the 
untracked state to the tracked one), which is OK 
and will result in zero chance of a spin flip when 
the positive random number is chosen. If no spin 
flip is to occur, skip to step [7] 

6. If a spin flip is to occur, compute the difference 
in potential energy between the currently tracked 
state and the other state, that is, by how much 
would an atom’s potential energy change if it were 
to flip from the tracked state to the untracked one 
at this point in space? 

A-Epot = b'lnitracked (A ^tracked (z) (2?) 

= QFijYln untracked U7. n tracked)MB |U(z) |. (28) 

Compare with the atom’s kinetic energy Ak; n . If 
Akin + AEp 0 t < 0, the spin flip is not (classically) 
energetically allowed. Continue tracking the origi¬ 
nal state and skip to step [7] 

Otherwise, modify the magnitude of the atom’s ve¬ 
locity v z , but not the direction, so as to conserve 
energy: 


Akin —t Akin + AAp 0 t (29) 


2 (Akin + AA pot ) 


(30) 


Finally, actually perform the spin flip—flipping en¬ 
tails simply noting that the other state is now being 
tracked, and does not involve modifying the wave- 
function in any way. 

7. Exponentially decay for one timestep the state not 
being tracked: 

lx) -t f 1 - ^Auntracked j |x) > (31) 


8. Normalize the wavefunction: 


lx) -» 


lx) 

lx)’ 


and proceed to step [I] 


(32) 


RESULTS 

In this section we compare our method to two oth¬ 
ers: the full spatial time-dependent Schrodinger equa- 
tiorj^] (TDSE), and the standard semiclassical method 
(which we refer to here as the Ehrenfest method) men¬ 
tioned in the introduction, whereby we calculate an ex¬ 
pectation value for the force on the atom using the in¬ 
stantaneous spin populations of the atom. In this section 
‘Monte Carlo wavefunction method’ or MCWF refers to 
our method. 

We first test our model on a simplified form of the 
Stern-Gerlach experiment, and subsequently test on a 
modified form of the Majorana problem which describes a 
ID magnetic trap with the possibility of Majorana losses 
near the center. 

The results shown for both the MCWF and Ehrenfest 
simulations are ensemble averages for simulations done 
with 10,000 atoms. The simulations begin with a station¬ 
ary, 20 pK 8 7 Rb atom that has a minimum uncertainty 
wavepacket: the spatial and momentum distributions are 
Gaussian with widths Ath and h/2X t h respectively, where 
A t h = h/\Z2 tuvEqT is the thermal de Broglie wavelength. 

Throughout this section we treat 87 Rb as a spin-^ 
atom with a Lande ^-factor gp equal to twice that of 
87 Rb’s actual spin-1 hyperfine groundstate. 

Population densities (and populations) for the TDSE 
are computed simply as the populations (and integrals) of 
its two internal spin states, in the basis of the local mag¬ 
netic field. The same quantities for the MCWF method, 
however, are computed not by reference to the internal 
state of the atom, but instead by counting what propor¬ 
tion of the atoms are being tracked as spin up or spin 
down. 


4 Simulated using the XMDS m differential equation solver. 
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FIG. 2. Results of the MCWF, TDSE and Ehrenfest meth¬ 
ods applied to a simple 1-dimensional Stern-Gerlach simu¬ 
lation. The shading shows the results of the TDSE simula¬ 
tion, the colored lines show the trajectories predicted by the 
MCWF method, and the black lines illustrate the trajectories 
predicted by the Ehrenfest method. Solid lines represent the 
mean trajectories, and dashed lines show one standard devi¬ 
ation from the mean. The plot at the top of the figure illus¬ 
trates the probability densities predicted by all three methods 
at the end of the simulation (t = 380 ps). Here solid colored 
lines are the results of the MCWF simulation, dashed colored 
lines represent the predictions of the TDSE simulation, and 
the black line is the prediction of the Ehrenfest method. 


The Stern-Gerlach experiment 

As discussed in the introduction, a prototypical exam¬ 
ple of spin-position entanglement is the Stern-Gerlach 
experiment. The original Stern-Gerlach experiment saw 
spin-half atoms passed through an inhomogeneous mag¬ 
netic field and after some propagation two distinct atom 
clouds were seen, implying the passage of two distinct 
trajectories, corresponding to the two distinct spin pro¬ 
jections along the direction of the field. Our method 
captures the essential physics of the Stern-Gerlach ex¬ 
periment in one dimension by simulating a spin-half atom 
beginning at rest in an equal superposition within a mag¬ 
netic field gradient along 2 . Figure [2] shows the results 
predicted by our MCWF method contrasted with TDSE 
and Ehrenfest solutions. 

Our initial wavepacket is centered at z = 7.8 pm in the 
magnetic field B (z) = (0,0 ,B' z z), with B’ z = 2.5Tm _1 . 
The wavepacket begins at rest, in a superposition of spin 
projection states (c-j- = cj, = 1/^2), which accelerate 
apart in the magnetic field gradient. This example shows 


perhaps the simplest case in which using an expectation 
value for the force on an atom is completely inapplicable. 
We see that the Ehrenfest approach results in the atoms 
remaining stationary, as one would expect given that the 
average of the forces the two components experience is 
zero. The shaded density plot shows qualitative agree¬ 
ment between the MCWF and TDSE simulations, with 
the respective spin components tracking similar trajec¬ 
tories. The inset shows the quantitative accuracy of the 
MCWF method, displaying good agreement between the 
density profile it predicts and that of the TDSE method 
at the final time step of the simulation. The noise ev¬ 
ident in both the Ehrenfest and MCWF simulations is 
statistical, and can be reduced by the addition of more 
simulation atoms. 


The Majorana problem 

Earlier we introduced the concept of a Majorana tran¬ 
sition and discussed the detrimental effect it can have 
during an ultracold atom experiment. In other work 
of ours, we want to model this effect using a full three 
dimensional gas simulation, modeling classical positions 
and momenta in order to scale to large numbers of atoms. 

Majorana derived a result [12] predicting the probabil¬ 
ity of a transition based on the ratio between the rate 
of change of the magnetic field direction and the Larrnor 
frequency. However this derivation was made for a sta¬ 
tionary atom in a dynamic field, which- unlike a moving 
atom in a spatially varying field does not develop spin- 
position entanglement. 

Moving into the center of mass frame of a moving 
atom in a spatially varying field is not sufficient to make 
the situation equivalent to Majorana’s assumptions, since 
he assumed a time-varying magnetic field only, with no 
spatial gradients that could give rise to forces. It was 
also derived in the asymptotic limit of a large longitu¬ 
dinal field—initially aligned with spin—which is steadily 
inverted in the presence of a constant transverse field. 
Nonetheless we consider a similar situation with a mov¬ 
ing atom subject to forces in a magnetic field gradient, 
and compare the spin flip probabilities it predicts to those 
predicted by Majorana for the effectively time-dependent 
magnetic field the moving atom sees at the field mini¬ 
mum. 

We apply the MCWF method to a simulation of an 
atom moving in a magnetic field gradient along a single 
spatial dimension, the results of which are shown in Fig¬ 
ures [3] and [4] The wavepacket is initially centered around 
2 = —50 pm in the magnetic field B(z) = (B x , 0, Biz), 
with B x = 105 nT and B' z = 2.5 Tin -1 . In figure [3[we 
have again plotted the trajectories predicted by each 
method and provided an inset displaying the probabil¬ 
ity densities at the final time. 

Figure [4] shows the evolution of the population in each 
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FIG. 3. The modified Majorana problem of an atom accel¬ 
erating in a 1-dimensional magnetic field gradient, simulated 
using the three methods: MCWF, TDSE and Ehrenfest. The 
density plot at the top of the figure is for t = 3.0 ms (as shown 
by the dashed line). The legend is as described in figure [5] 



FIG. 4. State population of the modified Majorana prob¬ 
lem in figure [3] here extended to t = 5.5 ms. The Ehrenfest 
approximation correctly predicts populations for the first min¬ 
imum crossing but subsequently fails for successive crossings. 


spin state over time. Here we show more simulation 
time to emphasize the failure of the Ehrenfest method. 
Again we have compared the results of the three differ¬ 
ent methods; the TDSE, the MCWF and the Ehrenfest 
methods. As before we note good agreement between 
the TDSE and the MCWF methods in the density plots, 
and observe the accuracy of the MCWF method in the 
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FIG. 5. The modified Majorana simulations over a range of 
transverse magnetic fields tests the technique over a range of 
spin transition probabilities. 


density profiles plotted in the inset. While the Ehren¬ 
fest trajectory initially agrees closely with the other two 
simulations, it diverges after the wavepacket has passed 
through the field minimum. The Ehrenfest populations 
agree well through the first crossing of the minimum and 
it isn’t until the trapped atom encounters another mag¬ 
netic field minimum that we really observe the failure of 
the Ehrenfest method in simulating state populations. 
This failure may be entirely attributed to the Ehren¬ 
fest method’s manifest failure to correctly simulate the 
wavepacket trajectories—most of its atoms are on an es¬ 
cape trajectory and will not see a second crossing of the 
magnetic field minimum. The stochastic nature of the 
MCWF method is evident in trajectory of the flipped 
state near the magnetic field minimum (figure [3| : here 
the population in the flipped state is small and so has 
large statistical noise. This noise decreases as the the 
flipped state population increases towards its asymptotic 
value. 

To ensure robustness of the MCWF method over a 
practical range of spin flip probabilities, in figure [5] we 
have plotted the predictions of the three methods over 
a range of transverse magnetic fields. As the transverse 
field strength increases the rate of change of the mag¬ 
netic field decreases and we expect to see a decrease in 
the probability of spin flip transitions. Each point in 
figure [5] is the final result of a simulation of the type 
shown in figure [3j These simulations all begin with a 
static 20 pK wavepacket of 10,000 s 'Rb spin up atoms 
at z = —29.8 pm. This is far enough away from the 
field minimum to be considered asymptotic. The atoms 
are allowed to accelerate through the magnetic field until 
they reach the antipodal point of the trapping potential. 
As expected we see the MCWF method agrees within its 
statistical uncertainty with the TDSE method. As be¬ 
fore the Ehrenfest method produces good results, but if 
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we were to allow the simulation to run longer we would 
soon see the method disagreeing with the MCWF and 
TDSE methods due to differing trajectories. 

DISCUSSION 
Computational costs 

The full many-body wavefunction is entirely infeasible 
to simulate for any more than a small number of atoms. 
Making the assumption that the atoms are in a product 
state allows on to simulate a single spatial Schrodinger 
equation per atom. Compared to this, our method has 
the computational cost of one ODE per atom instead of 
one PDE, making it tractable for large numbers of atoms. 

The results in the previous section do not demon¬ 
strate any such computational savings, since thousands 
of atoms were compared to a single TDSE run to demon¬ 
strate correctness. If one is content with only a single 
statistically representative outcome however, as is likely 
satisfactory for large thermal clouds of atoms, one need 
only perform one MCWF run. As in all stochastic meth¬ 
ods, repeated runs can be used to establish the variability 
of results at the cost of more computing resources. 

Applicability 

Our method of determining spin flip probabilities from 
population transfer allows us to more accurately model 
the classical trajectories of the simulated atoms than if 
we waited for population to be discarded, as in existing 
MCWF methods. We stop short of claiming, however, 
that computing spin flip probabilities in this way is at all 
correct for simulating other types of open systems with 
quantum jump methods. 

Indeed, it seems that if one were to try to apply this 
method to the common case of spontaneous emission of 
a photon from an atom, it would be the photon number 
states, not the atomic states, between which population 
transfer would have to be monitored. There is an extra 
layer present in this case compared to ours: the external 
environment is measuring photon number states, which 
in turn measure internal states of the atom. The time 
photon number states take to decohere, conditioned on 
photon detection with a large environment, is so short 
that even if one decided to use our method here, the 
post-jump exponential decay would be practically instan¬ 
taneous and indistinguishable from the immediate jumps 
these methods use. 

Another difference between our system and a decay¬ 
ing atom is that the positional states corresponding to 
each spin state are non-orthogonal (being mostly over¬ 
lapping Gaussians), unlike photon number states. This 
means that the assumption of constant, strong position 


measurements does not produce complete collapse of the 
spinor wavefunction, as any given position measurement 
provides only weak spin information. It is therefore 
natural that our model does not produce discontinuous 
changes in the system wavefunction. 


Higher spins and higher spatial dimensions 

Generalizing our method to higher spins entails defin¬ 
ing classical probability flows between each pair of states 
in a way that is consistent with the actual changes in pop¬ 
ulation at each timestep. Inferring classical probability 
flows is trivial for a two state system, but more difficult 
for higher spins due to the inevitability of interference. 

A recipe for inferring classical probability flows in a 
quantum system comprises a hidden-variables theory , de¬ 
fined by Aaronson as “a way to convert a unitary matrix 
that maps one quantum state to another into a stochastic 
matrix that maps the initial probability distribution to 
the final one in some fixed basis” [3T(. According to this 
definition, our method is a hidden-variables theory, with 
the spin projection quantum number m n being a hidden 
variable. Approaches to hidden-variables theories for dis¬ 
crete systems of more than two states exist, including the 
‘flow’ and ‘Schrodinger’ theories [5Tj . 

The only difficulty in generalizing to greater than one 
spatial dimension is knowing in which direction to mod¬ 
ify atomic velocities when kinetic energy is lost or gained 
during a spin flip. One solution might be to calculate the 
instantaneous force on a moving classical dipole, and to 
apply a velocity difference to the atom in the direction 
of this force. Another approach is to simply continue 
to apply velocity changes in the direction of motion in 
three dimensions. If the system is chaotic and there is no 
systematic bias making resulting trajectories statistically 
distinguishable from the correct ones, then such unphysi¬ 
cal velocity adjustments may nonetheless be satisfactory. 


CONCLUSION 

We have presented a one dimensional model for semi- 
classical simulations of spin-^ atoms in magnetic fields. 
Our model reproduces the salient features of the under¬ 
lying exact model, namely correct spin flip probabilities 
and statistically representative trajectories, which the 
Ehrenfest method fails to do. 

Our method is a great improvement over using Ehren¬ 
fest’s theorem for force calculation in semiclassical sim¬ 
ulations in contexts where superpositions of states sub¬ 
ject to different forces are non-negligible and long-lived. 
We intend to apply the method to simulations of forced 
evaporative cooling, in the hopes of further closing the 
gap between theory and experiment on the subject. 
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SUPPLEMENTARY INFORMATION: 


DERIVATION OF SPIN DECOHERENCE TIME r 


Inner product of two Gaussian wavepackets accelerating apart 


We want to know what the inner product of two equal width Gaussian wavepackets ipi and ij)j is as a function of 
time as they move apart with constant acceleration with magnitude : 




7*^ (x X r pl ) , - i 

e~4?e- dx ? 


(1) 


where 


and 


^rel(^) — 2 O'ijt 

jjl 

t) — —-Clijt 


h 

are the wavepackets’ relative position and wavenumber due to acceleration for a time t starting from rest, and: 


( 2 ) 

(3) 


C~' = 


e 2 ^ da; 


(4) 


is a normalization constant. Note that this expression holds for any number of dimensions — relative motion is only 
along one axis so the integrals in all other directions equal one. 

Let’s evaluate this integral by expanding the whole exponent into a polynomial in x, and completing the square in 
the exponent: 

4<t 2 
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da; 

da;. 


Now we’ve got a polynomial in x in the exponent and can complete the square: 
p(x) = x 2 - [x re \ - 2ia 2 k re i\ X + -x 2 el 

= X [#rel 2z<7 &rel] X — [#rel 2 i(T /Crel] — [#rel 2 id k re ij ~b ^^rel 


— 2 [^rel 2 i(T /c re l] ^ ^ [^rel 2i(J k re \j 
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= exp 
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( 12 ) 


(13) 


where c is a complex number. We recognize the remaining integral as a Gaussian integral with complex offset - it 
equals C~ x as in |4]). So it cancels C and we’re left with: 


{ipi{t)\il)j(t)) = exp 


gfj2 [ a ' re l 2*a fc re l] / |^ r 2‘ Z ' re l 


4a 2 ‘ 


(14) 
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FIG. 1. The computed decoherence time between the uif = ±1 states of the F = 1 hyperfine groundstate of 87 Rb in a 
125 Gem -1 magnetic field gradient. At small wavepacket sizes, pos ition separation dominates the decoherence (equation (211), 
whereas at larger wavepacket sizes velocity separation (equation ( |23| )) dominates. The simple expression (251 gives a good 
approximation to more complicated exact solution given by pOl 


over the domain, although the latter was not computed over the entire domain of the plot due to numerical overflow. 


8cr 2 


Expanding: 

('i/) i (t)\'ip j {t)) = exp 
gives us our result in terms of x re i and fc re i: 

= exp 

=> \(Mt)\ VbO))] = exp 


[•k’rel 4*tX ■Z'rel^'rel 4tJ /c rel ] ^ 2 *reli 
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^ 2 a h 2 
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(15) 


(16) 

(17) 


We see that for small wavepacket sizes <r, position separation dominates the decrease in wavepacket overlap, whereas 
for large wavepackets, velocity separation dominates. If it is known which of two terms in the exponent above 
dominates for a particular problem, here would be a good time to make the approximation of neglecting the smaller 
term. This will result in a simpler expression for the separation time below. 

Let’s substitute back in the expressions for x le \ and /c re i to get a polynomial in t in the exponent: 


|(V> i (t)|V’j(*))| = exp 
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(18) 


Time of separation 


Let’s treat this as an unnormalized probability distribution for the wavepackets still being overlapped, and thus 
define the average separation time t 13 = ( t) as: 


/o 00 *i(y , »(*)iy , j(*))i dt 
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This results in the following mouthful: 


7~ij exact — 


maaijKi [rf) 


( 20 ) 


where erfc is the complementary error function, A'i is a modified Bessel function of the second kind and 


77 = m 2 a 3 a 


is a dimensionless parameter. 


The above is not just inconvenient to calculate, in this form it is numerically difficult too, as the exponential 
overflows easily despite the overall expression being equal to a value that would not overflow. So let’s come up with 
an approximation. There are two regimes in which we can get very simple expressions for r, (7 -, they are in the limits 
of small and large wavepacket sizes. For small wavepackets (77 -C 1) we neglect the second term in ( [T?] ). This is 
equivalent to ignoring the velocity difference between the two wavepackets, and so gives us a time constant for the 


decrease in overlap purely due to separation in position space. Solving the integral (19) with this approximation gives: 


(27t 1 2 )4 1 a 

TuPOS ~ 2T(| )\ja ij 

(21) 

« 1.163. P?-. 

(22) 

V °*? 


Neglecting instead the first term in gives us the time constant for velocity separation, valid when 77 3 > 1 : 



h 

« 0.798-. 

mcraij 


(23) 

(24) 


Unless 77 ~ 1, one of these time constants will be much smaller than the other, and so the larger can be neglected if 
this is known to be the case. However, for an approximate expression that is valid for all wavepacket sizes, one can 
add the cubes of their reciprocal^] giving: 
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(25) 


A plot showing the accuracy of this expression compared to the exact solution is shown in figure [T] 

Note that for large wavepacket sizes, the single mode approximation may not be valid: a magnetic field with features 
smaller than the wavepacket size may induce spin transitions in some parts of the wavepacket but not others. The very 
short decoherence times shown in figure [l] for large wavepackets would therefore be unphysical. If the spin decoherence 
time is short compared to the timescale of spin transitions, this may unphysically suppress spin transitions via the 
quantum Zeno effect [32]. Caution should therefore be exercised not to choose a wavepacket size larger than the range 
over which the single mode approximation can reasonably be expected to hold. 

The simulations in our paper use the thermal wavelength at a particular temperature for the wavepacket size, which 
is the size of a wavepacket in a thermal gas shortly after a collision—considerably shorter than the time averaged 
wavepacket size, which is on the order of the mean free path [ 6 ]. 


1 If these time constants were for pure exponential decay, we could 

add their reciprocals directly. If they were for Gaussian decay, 

we could add the squares of their reciprocals. However one of our 


time constants is for Gaussian decay and the other is for decay 
which is quartic in t in the exponent, as per (jTsJ, so adding the 
cubes of their reciprocals (3 being close to the geometric mean 
of 2 and 4), gives a good approximation. 













